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In order to investigate the importance of pre-equilibrium dynamics on relativistic heavy-ion 
collision observables, we match a highly non-equilibrium early evolution stage, modeled by free- 
streaming partons generated from the Monte Carlo Kharzeev-Levin-Nardi (MC-KLN) and Monte 
Carlo Glauber (MC-Glb) models, to a locally approximately thermalized later evolution stage de¬ 
scribed by viscous hydrodynamics, and study the dependence of final hadronic transverse momentum 
distributions, in particular their underlying radial and anisotropic flows, on the switching time be¬ 
tween these stages. Performing a 3-parameter fit of the measured values for the average transverse 
momenta {p±) for pions, kaons and protons as well as the elliptic and triangular flows of charged 
hadrons with the switching time Ts, the specific shear viscosity rj/s during the hydrodynamic 
stage, and the kinetic decoupling temperature Tdec as free parameters, we find that the preferred 
“thermalization” times Ts depend strongly on the model of the initial conditions. MG-KLN initial 
conditions require an earlier transition to hydrodynamic behavior (at Ts ~ 0.13 fm/c) , followed 
by hydrodynamic evolution with a larger specific shear viscosity p/s ^ 0.2, than MG-Glb initial 
conditions which prefer switching at a later time (r^ ~ 0.6 fm/c) followed by a less viscous hy¬ 
drodynamic evolution with p/s ^ 0.16. These new results including pre-equilibrium evolution are 
compared to fits without a pre-equilbrium stage where all dynamic evolution before the onset of 
hydrodynamic behavior is ignored. In each case, the quality of the dynamical descriptions for the 
optimized parameter sets, as well as the observables which show the strongest constraining power 
for the thermalization time, are discussed. 

PACS numbers: 25.75.-q, 25.75.Cj, 25.75.Ld, 24.10.Nz 


I. INTRODUCTION 

The collision of two fast moving heavy nuclei produces 
a new form of matter, quark-gluon plasma (QGP). It 
consists of deconfined quarks and gluons, and in the 
early years QGP was thought to behave like a weakly- 
interacting gas. However, experiments conducted at 
the Relativistic Heavy Ion Collider (RHIC) and Large 
Hadron Collider (LHC) show strong collective flow, re¬ 
vealing that the QGP is in fact a strongly coupled liquid, 
with very low viscosity. The evolution of this strongly 
coupled QGP is well described by hydrodynamics, which 
successfully reproduced and even predicted the exper¬ 
imentally observed transverse momentum spectra and 
flow anisotropies of emitted hadrons. 

The validity of hydrodynamics relies on the matter be¬ 
ing close to local thermal equilibrium. Within a purely 
hydrodynamic approach that ignores details about how 
the system approaches equilibrium, the observed large 
hadron momentum anisotropy can only be explained if 
thermalization happens fast and the hydrodynamic ex¬ 
pansion does not begin later than about Ifm/c after the 
two nuclei impact each other [T] . This finding calls for a 
mechanism to explain how the matter produced in the 
collision can thermalize so fast. Recently, significant 
work has been done on modeling the thermalization pro¬ 
cess during the pre-equilibrium stage While these 
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studies have not yet fully explained the rapid thermaliza¬ 
tion, some of them have indicated that pre-equilibrium 
evolution can have a non-negligible influence on the final 
observables. 

In this paper, we return to the question of how quickly 
thermalization must happen for a hydrodynamic ap¬ 
proach to provide a successful description of experimen¬ 
tal data, by studying the weakly interacting limit of pre- 
equilibrium dynamics. Freely streaming partonic quanta 
can be considered as the extreme limit of a weakly inter¬ 
acting system and the diametrically opposite to a fluid 
dynamical description which requires strong coupling. 
By coupling a free-streaming pre-equilibrium stage to 
hydrodynamics and varying the switching time one 
can smoothly interpolate between a very strongly coupled 
(small Ts, early transition to hydrodynamics) and weakly 
coupled (large , late transition to hydrodynamic behav¬ 
ior) early evolution stage. We are interested in finding 
the largest value that is compatible with phenomenol¬ 
ogy. Replacing the non-interacting free-streaming par- 
ton stage in our model by a model in which the matter 
constituents interact more realistically with modest in¬ 
teraction strength should allow for an earlier transition 
to hydrodynamic behavior. For this reason we expect our 
approach to yield a robust upper limit, for a fixed viscos¬ 
ity during the subsequent hydrodynamic stage, for the 
time at which the matter created in the collision must 
have reached a sufficient degree of local equilibrium to 
be described hydro dynamically. The dependence of this 
upper limit on the viscosity in the hydrodynamic stage 
will also be explored. 
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It was initially thought that substituting the earli¬ 
est stage in a hydrodynamic evolution model by a free- 
streaming gas would delay the buildup of thermodynamic 
pressure and thus reduce the finally observed collective 
flow [1]. This ignored, however, the fact that in a spa¬ 
tially inhomogeneous medium free-streaming (or, for that 
matter, any kind of pre-equilibrium evolution) gener¬ 
ates strong position-momentum correlations which, upon 
thermalization, lead to strong initial flow in the hydro- 
dynamic stage. We show here that this “pre-flow” actu¬ 
ally increases the finally observed radial flow, with conse¬ 
quences that are opposite to the expectations reported in 
Refs. [1]. Here we explore the evolution of the widely im¬ 
plemented MC-KLN [T0l[TT] and MC-Glauber [12] initial 
conditions at LHC energies, allowing them to be evolved 
by free-streaming for a time Tq before switching to a vis¬ 
cous hydrodynamic description for the rest of the evolu¬ 
tion. 

Our work focuses on the effects brought by free- 
streaming on both the hydrodynamic initial conditions 
and the final observables. We mainly focus on the MC- 
KLN model which provides a complete prediction for the 
initial gluon distribution, not only in space but also in 
momentum. However, we show that for massless partons 
moving with the speed of light the shape of the initial mo¬ 
mentum distribution is irrelevant as long as it is locally 
isotropic. This allows to apply our description also to 
MC-Glb initial conditions although that model makes no 
prediction per se about the initial parton momentum dis¬ 
tribution. Free-streaming evolves the initial conditions 
from an initial parton formation time tq (which is taken 
to be very close to zero) to the switching time when we 
switch to a near-equilibrium hydrodynamic description. 
The sudden transition to approximate local equilibrium 
is implemented by applying the Landau matching proce¬ 
dure. By tuning the switching time, we can enforce fast 
thermalization by setting Tq—To^ or slow thermalization 
by setting ^tq. The hydrodynamic initial conditions 
obtained from the Landau matching procedure vary with 
the switching time, enabling an investigation of the in¬ 
fluence of Ts on the final observables. The hydrodynamic 
evolution is performed with the code VISH2+1 El 13, 
without hadron cascade afterburner. For the hydrody¬ 
namic evolution, we use a constant specific shear vis¬ 
cosity Tj/s. Freeze-out is implemented at a fixed kinetic 
freeze-out temperature T^ec, followed by a Cooper-Frye 
procedure with full resonance decay cascade to convert 
the hydrodynamic output into final stable particle spec¬ 
tra. 

In Section [H] we outline the free-streaming evolution 
in the pre-equilibrium stage and describe the Landau 
matching procedure. Its consequences on the hydro- 
dynamic initial conditions are discussed in Section in 
Section [LV] shows how the hydro dynamical evolution re¬ 
sponds to initial conditions generated at different switch¬ 
ing times. A difficulty related to the conversion of par- 
tons to hadrons that arises from a late switching time 
is discussed and resolved in Section [Vl In Sections I VII and 


|vnl the energy flow anisotropy and hadron mean trans¬ 
verse momenta are constructed to illustrate how final ob¬ 
servables change with switching time. Finally in Section 
1^ we introduce a multidimensional parameter search 
procedure to systematically study the preferred ranges of 
Ts, Tj/s and Tdec- For both MC-KLN and MC-Glb ini¬ 
tial conditions, both with and without a free-streaming 
pre-equilibrium stage before the onset of hydrodynamic 
behavior, we determine the best-fit parameters and their 
uncertainty ranges, and discuss their relative quality of 
describing the data. Conclusions are presented in Section 
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II. FORMULATION OF FREE-STREAMING 
AND LANDAU MATCHING 

The evolution of partons in the free-streaming model 
is described by the collisionless Boltzmann equation 

p^d^f{x,p) =0. ( 1 ) 

We work in Milne coordinates and write f{x^p) = 
/(3^±,^), with longitudinal proper time 
r = space-time rapidity r]s = ^ \n[{t-^z) / (t—z)], 

and rapidity y = ^ln[{E^Pz)/{E—pz)]. We assume 
massless partons for which E=\p\ = The 

collisionless Boltzmann equation is easily solved 
analytically, relating the final parton distribution 
f{x±,r]s,rs;p±,y) to the /(x_l, tq;^) by a spatial 
coordinate shift, keeping the p± distribution unchanged. 
For massless partons one finds: 

f{x±, ris, Ts-,pj_,y) = fix±-{Ts-To)p±, 11s, To] p±,y), (2) 

where =p±/p± = (cos sin^p), with <pp being the 
azimuthal angle of p± in the plane transverse to the 
beam. The initial distribution /(cc^, 77 ^,tq; px,^) is as¬ 
sumed to be locally isotropic in transverse momentum, 
i.e. independent of /o = f{x±,r]s^ro;p±,y). 

To initialize the hydrodynamic we must decompose the 
energy momentum tensor in hydrodynamic form. In 
the free-streaming stage {x) can be obtained from the 
solution of Eq. <§ for the parton distribution as 

Tf^‘'{x±,ris,Ts) = j^^ j ^P^p‘'fix±,ris,Ts]P±,y),i3) 

where for massless particles p^ = E =\p\ and ^ is a degen¬ 
eracy factor. We assume longitudinal boost-invariance 
and restrict the dependence of f{x,p) on r]s and y as 
follows nmn]: 

f{x_i,ris,T] p±,y) = ~^ f{x_i,T] p±,y). (4) 

Again, for massless partons m± = = p ±• For 

a boost-invariant system it is sufficient to know in 
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the transverse plane at z = 0 (i.e. 77s=0) where 


7/5=0, r) = 


( 27 r)^ T 


2 |*c>o /*7r 

- dp± 

T Jo i-TT 


d(pp 


P^P''f{x±,T;p±,y) .( 5 ) 


y=0 


The factor 1 /r is characteristic for systems with boost- 
invariant longitudinal expansion. 

Inspection of this formula shows that for massless par- 
tons the evolution of ^s=0, r) does not depend 

on the underlying parton momentum distribution. For 
massless partons that are initially distributed locally 
isotropically in p± the spatial distribution of the energy 
momentum tensor at time r depends only on 

its initial spatial distribution at time tq, but not on the 
px-distribution (which may or may not depend on x±). 
One sees this by observing that the integral on the right 
hand side of Eq. <§ can be rewritten as 


pOO ^TT 

/ dp± / d4>p p^p''f{x±,T;p±,y) 

Jo J-n 

pOO pTT 

/ Pldp± d4>pp'^p‘'f{x_i-{T-TQ)p±,TQ;p±,0){6) 
Jo J- 7 T 


where ^ 


depends only on Hence 

y=0 


T^'"(®x,7?s=0,t) = t f d(pp p>"p''F{x±,T;(pp), (7) 

J — IT 

where 


F{x^, t; (j)p) = Fo{x±-{t-to)pj_) 

P±dp±f{x±-{T-To)p±,To;p±,0) ( 8 ) 

is independent of how /o depends on the magnitude of 
p_L, and Fo{x±) denotes the spatial distribution function 
(integrated over momenta) at r = ro. 

At the switching time r^, the solution 0 for is 
decomposed in viscous hydrodynamic form: 



= eu^u^ -{VF n)A^^ + tt^A (9) 


Here e and V are the energy density and pressure in 
the local fluid rest frame (LRF), H is the local bulk 
viscous pressure, is the shear pressure tensor, and 
is the local fluid velocity. The projection operator 
projects on the spatial coordinates in 
the LRF, and the spacetime metric in Milne coordi¬ 
nates is given by g^^ = diag(l, —1, —1, —1/r^). 

The Landau matching condition defines the fluid rest 
frame velocity as the time-like eigenvector of and 
the energy density e as its eigenvalue: 


( 10 ) 

with = 1 . If such a solution exists m it is unique 

since has at most one time-like eigenvector. 

Landau matching conserves the system’s total energy, 
but not its entropy. Due to the absence of collisions. 


Eq. 0 conserves entropy. The solution 0 therefore 
yields the same entropy at and tq. After matching, 
however, the entropy density s = dVjdT is related to 
the energy density e and pressure V by the thermalized 
fluid’s equation of state (EOS) V = V{e). This implies 
that in general the total entropy S of the system in¬ 
creases discontinuously at the switching time Tq. This 
sudden increase is the consequence of the assumed sud¬ 
den thermalization of the system that is implicit in the 
Landau matching procedure. For successful phenomenol¬ 
ogy, we have to normalize the entropy density profile after 
Landau matching such that, upon completion of the dy¬ 
namical evolution, it correctly reproduces the observed 
final multiplicity dNch/dy. This is done for central col¬ 
lisions, and the predicted impact parameter dependence 
of the final dNch/dy is then taken as an argument for or 
against the validity of the initial model used to generate 
the initial conditions. In spite of the entropy jump at 
Tg, the normalization of the entropy density profile after 
Landau matching has a one-to-one relation with the nor¬ 
malization of the initial distribution function. Since the 
entropy jump depends on the chosen value of the switch¬ 
ing time, this initial normalization also depends on r^: 
to preserve the same dNch/dy at the end of the hydro- 
dynamic evolution thus requires a renormalization of the 
initial distribution function when is varied. 

After finding the energy density, the thermodynamic 
pressure is given by the EOS 7 ^ = T’(e) of the thermalized 
liquid. The dynamically induced bulk viscous pressure is 
then reconstructed from by using the identity 

n = -iiv(A^,rn- 7 >. (11) 

Finally, the shear pressure tensor is obtained by con¬ 
tracting with the double projection operator A^^ = 

1 (a;(A^+a^a-) - |aa‘-a«^, 

( 12 ) 

Alternatively, 


^ ^ ^)A^^ ( 13 ) 

where is from Eq. 0’. 77^ and e from Eq. ( [Tq| ), V 
from EOS and H from Eq. 

The initial conditions for e, T’, H and obtained 
from Eqs. (10)-(12) are more realistic than those used in 
many recent implementations of hydrodynamics where 
any evolution of the system of the system between tq 
and the hydrodynamic starting time Tq is ignored and 
the hydrodynamic stage is initialized with zero trans¬ 
verse flow and an ad hoc guess for the viscous pressure 
components. We will see that the non-zero flow veloci¬ 
ties and well-defined non-zero initial conditions for H and 
resulting from the Landau matching procedure ( p!Q| ) 
have important consequences for the subsequent hydro- 
dynamic evolution and final observables. In this work, 
we focus on the effects of shear viscosity on the evolu¬ 
tion of the collision systems. We implement our assump¬ 
tion of vanishing bulk viscosity C by setting C/s = 10“^ 
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FIG. 1. (Color online) Contour plots of the local energy density multiplied by Ts in the transverse plane, for a single fluctuating 
Pb+Pb collision event of 10% - 20% centrality at ^/s = 2.76 ATeV. To show the edge of the hreball more clearly, white is used 
for Tge below 10“^ GeV/fm^, and dark red is used for r^e > lOGeV/fm^. 


and then using the second order evolution equation m 
{D = uf^d^) 


nn = --{n + (:e)-ln^dJ^uA (i 4 ) 

to evolve the bulk viscous pressure 11 dynamically to zero, 
on a microscopic time scale given by the bulk relaxation 
time rn = ^ [IZ]. 

III. HYDRODYNAMIC INITIAL CONDITIONS 
AFTER A FREE-STREAMING 
PRE-EQUILIBRIUM STAGE 

In this and the following sections where we investigate 
the qualitative effects of a free-streaming pre-equilibrium 
stage on the hydro dynamical evolution and final observ¬ 
ables, we focus on MC-KLN initial conditions mm- 
We will return to the MC-Glauber model in Sec. Eul 
Compared to hydrodynamics, free-streaming dilutes 
the local energy density much faster. Because there are 
no collisions, signals carried by the massless partons move 
with the speed of light instead of the smaller drift veloc¬ 
ity that would characterize an interacting medium (for a 
thermalized medium this would be the speed of sound). 
In Fig. the local energy density just after switching to 
hydrodynamics is shown at different switching time r^. 
In the first panel, Landau matching is implemented at 


the matter formation time tq (taken as ro=0.01 fm/c). 
At this time, the matter distribution features many hot 
spots in the transverse plane, reflecting the fluctuating 
nucleon positions that, through the nuclear thickness 
function, control the saturation momentum and thus the 
density of the produced gluons in the MC-KLN model. 
As the switching time increases, the bumps in the en¬ 
ergy density spread and gradually dissolve as a result of 
the free-streaming of the partons. The initially bumpy 
energy density profile becomes much smoother and less 
deformed, resulting in decreasing spatial eccentricity co¬ 
efficients Sn- 

The spatial eccentricity coefficient of harmonic order 
n at the beginning of hydro dynamical stage is usually 
defined as iniisn] 


£n — ' 


An^ri 


f d^r± e(r±) 

f (Pr± e{r±) 


n > 1 , ( 15 ) 


where e{rj_) is the LRF energy density obtained from 
Eq. and the minus sign ensures that the angle 
points to the direction where energy density falls fastest. 
(For fluctuating events, the energy density profile must 
be re-centered to the origin in the transverse plane before 
calculating the eccentricities.) However, if the initial con¬ 
ditions feature a non-zero initial hydrodynamic flow pro¬ 
file, as is the case in Eq. ( [Tq| ), flow anisotropies cause the 
Lorentz contraction factor 7 between the local and global 
rest frame to depend on the azimuthal angle (j). In the 
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laboratory frame, the initial energy density is thus better 
characterized by eccentricity coefficients calculated with 
a modified prescription using the energy density in the 
laboratory frame: 

£n{Ts) = (16) 

d^a^(x) T^^(x) Uiy(x) r'2^ 
f d^r±j(r±)e(r_L)r2e^"'^ 

/ d^r_L7(r_L) e(rj_) r" ’ ^ 


where (Pcr^ is the normal vector on the switching hyper¬ 
surface of constant Tg. Now points in the direction 
of the steepest descent in the lab frame. If need be, the 
modified eccentricity definition (16) can also be used for 
different switching surfaces. 


nite mean free path A^fp = oo), the Knudsen number 

Kn = (17) 

-^macro 

where Lmacro is the characteristic macroscopic length 
scale of the system, is infinite. This tells us that, even 
if the initial momentum distribution were thermal, the 
system would evolve further and further away from local 
thermal equilibrium. The inverse Reynolds number uses 
the hydrodynamic decomposition to describe how far 
away a system is from local thermal equilibrium. For a 
non-equilibrium pressure caused by shear viscosity, it is 
defined as the ratio between the scalar charac¬ 

terizing the magnitude of the shear stress and the thermal 
pressure V: 



FIG. 2. (Color online) Event-averaged spatial eccentricity as 
a function of switching time for an ensemble of 400 fluctuating 
Pb-hPb collision events of 10% - 20% centrality at yC = 
2.76 A TeV. 

As shown in Fig.[^ the event-averaged ellipticity 62 and 
triangularity 63 at the beginning of the hydrodynamic 
stage decrease as the length of the free-streaming period 
increases. However, for > 6 fm/c, they increase again. 
Figure [^illustrates why this happens: around 5 fm/c the 
free-streaming fireball begins to disintegrate and even¬ 
tually separate into multiple pieces. This disintegration 
happens due to the absence of interactions between the 
partons which would otherwise keep the fluid together. 
So the matter distribution becomes more eccentric for 
large switching times, and the eccentricity coefficients 
increase accordingly. 

Free-streaming also drives the system out of equilib¬ 
rium. For free-streaming (which corresponds to an infi- 


,_i _ 


= 


V 


(18) 


For an analytic estimate Eq. (18) is inconvenient since 


the thermal pressure V is related to the energy density 
e (whose initial profile can be calculated analytically in 
terms of the distribution function /) only through nu¬ 
merical lattice QCD calculations. This complication can 
be avoided by using a slightly different definition for the 
inverse Reynolds number: 


,_i _ 


R-^ = 




p + n ■ 


(19) 


For a conformal EOS, V = \e and 11 = 0, this defini- 
tion agrees with Eq. (18), but if H 7 ^ 0 it allows for the 


following analytic computation of its initial value. 

R~^ is a local parameter whose initial value for the 
case at Tq = tq can be calculated analytically, provided 
the initial momentum distribution is locally isotropic in 
the transverse plane and the system is boost-invariant. In 
this case, inserting the definition (§ into ( [T^ to obtain 




g d?p f g cRp' 


f gd-^p f 

J {2tt)^p° J 


(27r)y° J (27r)y^ 

gd^p f gd^p 




^ f gd-^p f 

J {2ttYp^ J 


{2tt)^P 


,/o 


2^/2 


ip-p') -op p‘ 


f{p)f{p'), 

( 20 ) 


using local transverse momentum isotropy to recast 
Eq. @ for massless partons in the transverse plane at 
r]s = 0 into 


f{p)\ 


To,r]s=0 


sjy) 

ToP±' 


fip±), 


( 21 ) 


rewriting f = / dyd^p±, and using the identity 
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FIG. 3. (Color online) Contour plots for the inverse Reynolds number R~^ (19) in the transverse plane just after the system is 
switched from free-streaming to hydrodynamics, for several choice of switching time for a single fluctuating Pb+Pb collision 
event of 10% - 20% centrality at ^/s = 2.76 ^TeV. See text for discussion. 


for massless on-shell gluons we find 




lro,r7s=' 




dyd^p^dy'd^p'^ 

(27r)6ro2p_Lp'_L 


1 2 /2 


{P P ) -oPP 


fiPx)fip'±)Siy)d{y') 


^ f dpj_dp'j 

To J (27r)6 

j d4)pd4>'p 


X 


^p\p'lf{pi-)f{p'x 


sin (ppsm 


y [{sin Op cos (l)p sin Op cos (l)p (22) 


+ sin Op sin (pp sin Op sin pp + cos Op cos O'pY — _ ,_q- 


In the second line, p was decomposed as p = 
p(sin Op cos pp^ sin Op sin cos Op) using spherical coordi¬ 
nates, with the polar angle Op {Op) being related to the 
pseudorapidity r] {r]') (which for massless partons agrees 
with their rapidity y {y')) through 


where C = ^ / j:^p\dp^f{pA_)- 

The value of 7^+11 (the denominator of R~^) at r = ro 
and = 0 can be found from Eq. 0 (all quantities 
evaluated at r = tq and 77 ^ = 0 ): 


V 


•n = — 

3 ^ 


= — 
3 


/ 


3(27r)3To 


/ 


gd^p , . f{p) 

27r 


(27r)^p* 

d?pi 


pif{pn = ^C. (25) 

p±_ d 


Combining Eqs. ( |24[ ) and (25), the inverse Reynolds num¬ 
ber R~^ (19) at the matter formation time tq is seen to 
be independent of position in the transverse plane and 
equal to 


R-^ 


TO ,^ s =0 



(26) 


COS Op = tanhr^, sin Op = 1/ cosh 77 , (23) 

since 77 = — In [tan {Op/2)]. We see that for massless parti¬ 
cles the angular integral can be factored from the integra¬ 
tion over p±. At y = y' = 0^ Op = 0p = and the integra¬ 
tion over azimuthal angles pp and p'p is easily performed, 
giving the result 27r^/3. Thus 

277 ^ 

(24) 


This non-zero value for the initial inverse Reynolds num¬ 
ber is caused by the anisotropy in the initial pressure 
whose longitudinal component vanishes due to the as¬ 
sumed zero width of the initial rapidity distribution 
6{y — r]s)‘ It is consistent with the value obtained in 
anisotropic hydrodynamics which gives for a conformal 
system [ 22 ] : 

(V^f^^)^hydro ^ /3 1-Vl/Vx , . 

V V21 + Pl/(2Px)’ ^ ’ 
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and thus agrees with our results at 7^^ = 0. 


The value of R~^ at the beginning of the hydrodynamic 
stage increases if the matter is allowed to free-stream for 
a finite time before switching to hydrodynamics. This 
is caused by the appearance of additional anisotropies 
in the transverse plane. As the partons free-stream and 
the matter expands outward, the transverse momentum 
distibution becomes locally anisotropic, especially in the 
regions near the outer edge of the fireball where the mo¬ 
menta point predominatly outward. This breaks the fac¬ 
torization in Eq. (23) of the angular integral from the 
one over the magnitude of p± , and the transverse profile 
of R~^ can no longer be calculated analytically. Fig. 
shows how the inverse Reynolds number R~^ just af¬ 
ter switching varies with r^. The constant value of R~^ 
shown in the first panel of this figure (where we switch 
directly at tq) reflects the result ( [^ . As increases 
in the following panels, R~^ first becomes large near the 
edge, but later also in the core of the profile, indicating 
that the entrie system is moving farther away from local 
thermal equilibrium if it is allowed to free-stream longer. 


One particular feature of the hydrodynamic initial 
conditions obtained by Landau matching after free- 
streaming is rather strong initial flow. Its radial part 
can be quantified by the mean radial velocity 



T, (fm/c) 


FIG. 4. Event-averaged initial mean radial velocity as a 
function of switching time r^, for an ensemble of 400 fluc¬ 
tuating Pb-hPb collision events of 10% - 20% centrality at 
2.76ATeV. 


IV. EFFECTS OF FREE-STREAMING ON 
HYDRODYNAMICAL EVOLUTION 


W} = 


/ (fr±-f{r±)v±{r±)e{r±) 
f d^r±'y(r±) e(r±) 


(28) 


where the curly bracket stands for the single-event av¬ 
erage over the transverse plane. The e and v± profiles 
are taken from the Landau matching results. For event- 
by-event hydro dynamical runs, the event-averaged initial 
radial velocity (v±} is obtained by averaging the value 
(28) just after Landau matching over the events: 




1 




events 

E 

i=l 


(29) 


Fig-i shows this quantity as a function of switching 
time. It initially rises very quickly, reaching 25% of the 
speed of light already after 1 fm/c, and continues to grow 
at an approximate rate {aj_) ~ = 0.13 c^/fm over 

the next 5 fm/c. We reiterate that we expect a free- 
streaming pre-hydro dynamic stage to yield the largest 
possible pre-equilibrium effects on the subsequent hydro- 
dynamic evolution. These will be discussed in the fol¬ 
lowing sections; some of the trends we will observe in 
these sections may manifest themselves in weakened form 
when free-streaming will be replaced by more realistic 
pre-equilibrium dynamical models in future studies. 


A. Radial flow 


Pre-equilibrium evolution endows hydrodynamics with 
significant initial flow but a reduced initial spatial eccen¬ 
tricity. The interplay between these two tendencies con¬ 
trols the radial flow and its anisotropies that are finally 
observed in the measured hadron momentum distribu¬ 
tions. For hydrodynamic simulations starting at differ¬ 
ent switching times, the initial conditions are rescaled 
to guarantee constant final total energy per unit space- 
time rapidity, dE/dijs^ as calculated from a single-shot 
hydrodynamic run with an event-averaged initial profile 
without pre-equilibrium evolution that has been tuned 
to reproduce the measured total charged hadron multi¬ 
plicity density dNch/dr] in central Pb+Pb collisions at 
^/s = 2.76 A TeV. The final dE/drjs is calculated on the 
freeze-out surface as 



(30) 


dE/dr]s is roughly equal to the total multiplicity den¬ 
sity dN/dr] multiplied by the mean transverse mass for 
charged hadrons, Note that by varying initial 

parameters keeping dE/drjs fixed, we allow dNch/dr] to 
vary: if the parameter change leads to an increase in 
radial flow, increases and dNch/dr] will decrease 

accordingly. The reason for using in this work the fi¬ 
nal energy dE/dr]s rather than dNch/dr] for rescaling the 
initial distribution for each fluctuating event will be ex¬ 
plained in Sec. |V[ 









We expect the increase with of the initial average 
radial flow shown in Fig. to manifest itself in a flow- 
induced blue shift of the finally measured hadron p±- 
distributions, caused by the hydrodynamic radial flow 
on the freeze-out surface. For a single event, the final 
hydrodynamic radial flow can be characterized by the 
average radial velocity of the fluid cells on the freeze-out 
surface 


/sfo e 


(31) 


In the second equality we used that our freeze-out surface 
has constant temperature T^ec = 120 MeV and, therefore, 
constant energy density edec • The value vfo of this quan¬ 
tity obtained from Eq. (31) for a single hydrodynamic 
run with an ensemble-averaged MC-KLN initial profile 
for Pb+Pb collisions at y^ = 2.76 74TeV and 10% - 20% 
centrality is shown by the black squares in Fig. as a 
function of switching time r^. As expected from Fig. 
it increases with r^, contrary to what we observed ear¬ 
lier in simulations that ignored pre-equilibrium dynamics 
and started the hydrodynamic evolution with unevolved 
density profiles and zero transverse flow at the same Tq 
(shown as red circles in Fig. [^. 



FIG. 5. (Color online) Vfo as a function of switching time 
Ts , obtained from a single-shot hydrodynamic event with and 
without pre-equilibrium evolution, for the same ensemble- 
averaged MC-KLN initial profile. The solid (dashed) line cor¬ 
responds to including (excluding) pre-equilibrium flow before 
beginning of the hydrodynamic evolution stage at Ts . 


B. Anisotropic flow 

The total momentum anisotropy 

, _ / d^rj_{T^^ -Tyy) 
~ f d^r±(T^^ + Tyy) ’ 


( 32 ) 


is directly and monotonously related to the elliptic flow 
of all charged hadron integrated over p_L [24l |25] . This 
quantity includes the collective flow anisotropy generated 
by pre-equilibrium dynamics as well as a contribution 
from the anisotropy of the local momentum distribution 
reflected in [8]. The collective flow part of this mo¬ 
mentum anisotropy is captured by 


^ jd?r^{T-- + Tyyy 


where 7]^^ = eu^u^—VA^^ is the ideal fluid part of the 
enery momentum tensor. Starting hydrodynamics at dif¬ 
ferent switching times causes both quantities to saturate 
at different values. At large times, velocity shear effects 
die out and 11, become small 0, and hence Cp and 
e'p approach each other (see Fig. 1^. Smaller saturated 
values of and indicate weakei^nal anisotropic flow. 
Studying how the saturated values of these quantities 
change with switching time thus may help to constrain 
that parameter. 

In order to calculate the ensemble-averaged total mo¬ 
mentum anisotropy (e^), we first rotate the transverse 
components of in each event to maximize the mag¬ 
nitude of e'p and then sum over events. For (e^) we pro¬ 
ceed similarly with the ideal fluid part of Fig. 

shows how these ensemble-averaged (e^) and (e^ evolve 
during the hydrodynamic stage, for several choices of the 
switching time, (e^ always starts out with zero mag¬ 
nitude at r = Ts since the initial parton momentum dis¬ 
tribution at To is isotropic, and this isotropy of the spa¬ 
tially integrated momentum distribution is preserved by 
free-streaming. In contrast, starts out after Landau 
matching with non-zero magnitude, as seen in Fig. [^. 
This is due to anisotropies in the space-momentum cor¬ 
relations that were generated in the pre-equilibrium stage 
and that, after Landau matching, manifest themselves in 
anisotropies of the hydrodynamic flow profile u^{x). Af¬ 
ter the onset of hydrodynamic evolution, (e^) initially 
increases quickly, rapidly approaching e^, and then sat¬ 
urates. For small switching times < 2 fm/c, (e^ has 
enough time to fully develop before freeze-out, and the 
saturated values show little sensitivity to r^. However, 
for larger switching times Tq > 2 fm/c, (e^ saturates at 
lower values that decrease rapidly with increasing r^. If 
thermalization happens very late, hydrodynamics is no 
longer able to reach the same degree of final momentum 
anisotropy as obtained for early thermalization. This 
agrees with earlier findings in [U [23] . 


V. ISSUES OF PARTON-HADRON 
CONVERSION 


As shown in Sec. IIF free-streaming decreases the local 
energy density. As the free-streaming time Tq increases, 
more cells fall with their energy density below the de¬ 
coupling value Cdec even before the hydrodynamic stage 
starts. These partons will not be able to thermalize. The 
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FIG. 6. (Color online) Time evolution of the ensemble-averaged total hydrodynamic flow anisotropy (e^) (a) and total mo¬ 
mentum anisotropy (e^) (b) (see text for definitions) for different switching times, for an ensemble of 400 fluctuating Pb+Pb 
collision events of 10% - 20% centrality at ^/s = 2.76 ATeV. 


sketch shown in Fig. illustrates this phenomenon. In 
this figure, the thick colored horizontal lines represent ar¬ 
eas occupied by partons. Blue segments stand for regions 
with energy density below Cdec- After Landau match¬ 
ing, these regions would be outside the freeze-out sur¬ 
face Ufo that separates the thermalized fluid from free- 
streaming particles. The particles in those regions “de¬ 
couple” instantaneously on a surface T^out that is part 
of the T = Ts Landau-matching surface. Partons emitted 
from Sout neither thermalize nor evolve hydrodynami- 
cally. The red segment labeled as Hhydro indicates the 
region with e > edec- After Landau-matching, cells on 
Shydro bccomc part of the thermalized fluid, so Shydro 
forms the initial condition surface from which the hy- 



FIG. 7. (Color online) Illustration of would-be fluid cells in¬ 
side (red) and outside (blue) of the freeze-out surface dur¬ 
ing free-streaming in a (l+l)-dimensional space-time diagram 
where the horizontal axis represents the transverse plane. 


drodynamic evolution starts. After hydrodynamic evo¬ 
lution, these cells decouple into free-streaming particles 
once they reach the freeze-out surface Ufo, i.e. once their 
local energy density drops below Cdec- 

At r = To, only a small region near the edge of the fire¬ 
ball is on Hout- When free-streaming starts, the system 
expands both longitudinally and transversally, leading to 
a decrease of the local energy density, so T^out grows and 
Shydro shriuks. At very large switching times ^ 8 fm/c, 
Hout covers the entire r = Tq surface, and hydrodynamic 
evolution would never even start. 

The non-thermalized partons on Sout hadronize di¬ 
rectly from the free streaming stage. Contrary to the 
partons at Hhydro whose hadronization happens only after 
hydrodynamic evolution on Sfo and is described by the 
Cooper-Frye formula [26], partons on Sout do not have 
thermal distributions and we cannot use the Cooper-Frye 
formula to convert them to hadrons. Furthermore, most 
of the partons generated from the MC-KLN model are 
soft, with characteristic momenta <1—2GeV. No reli¬ 
able models exist to convert such soft partons to hadrons. 
However, even though we cannot follow their evolution 
into final hadrons of well-defined mass and flavor, we can 
still follow their momentum and energy. So in the follow¬ 
ing sections, we propose to use the system’s total energy 
flow distribution, instead of the momentum distributions 
of identified hadrons, to investigate radial and anisotropic 
flow for late switching times. We will trust our model’s 
predictions for final identified hadron flows only when the 
contributions from T^out can be neglected. 

Our inability to hadronize the contribution from Sout 
causes a problem for the rescaling of the initial profile. 
In common practice, the initial entropy profile is rescaled 
such that the calculated dNch/drj matches the experimen¬ 
tal measurement. This procedure is problematic if their 
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is a significant contribution of final hadrons emitted from 
Sout that cannot be included. Our way to work around 
this problem is to normalize the final energy on the trans¬ 
verse plane to a “standard” value. This “standard” final 
energy is obtained from a smooth fireball evolved with 
single-shot hydrodynamics starting at 0.6 fm/c without 
pre-equilibrium dynamics, with parameters matched to 
reproduce the experimental dNch/dr]. Using longitudinal 
boost-invariance, i.e. the fact that the integrand can only 
depend on the rapidity difference y—r]s^ the total energy 
per unit rapidity on the freeze-out surface is given by 


dE 

dy 


9i 


= E 

i 

= E (27r)3 


[ (fp± [ p'"(r(T^i{n^Pv)fi 

J J-Eto 


[ 


drjs 


dE 

dr]s 


(34) 


where = ( 1 , 0 , 0 , 0 ) is the temporal unit vector in the 
lab frame and the sum runs over all species. The parton 
distribution for hadron species i at the freeze-out surface 
can be written as 


fi — fo,i T ^/shear,i T ^/bulk,i: (3^) 

where /o,i is the local equilibrium distribution for hadron 
species i while (^/shear,i and (^/buik,i) are the shear and 
bulk viscous corrections to /o,i, accounting for the sys¬ 
tem’s deviation from local thermal equilibrium. For the 
shear correction we use the ansatz [211 EH] 

5fshear,i = 

For the bulk correction we use the following expression, 
derived from the 14-moment approximation for particles 
with Boltzmann statistics [29] : 

^/bulk,i — /Ojifl- T -^ 0 ,i T ] 5 

(37) 


For a non-interacting hadron resonance gas the coef¬ 
ficients Bo{T), Do{T) and Eo{T) were calculated in 
Ref. [30] in the Boltzmann limit. 

Although in our work here the bulk viscous pressure 
approaches zero on the short bulk relaxation time scale 
Tn, its initial value from the Landau matching is large, 
and the bulk viscous correction (5/buik,i remains signifi¬ 
cant over the early part of the hydrodynamic freeze-out 
surface sketched in Fig. To match the last two lines 
of Eq. (34) it is therefore important to include (^/buik,i in 
the definition of the distribution function 


^ For a quantitatively accurate matching one actually must cor¬ 
rect the ansatz ( |37| and the calculation m of its coefficients 
Bo(T), Do{T) and Eo(T) for quantum statistical effects. We 


Even if we do not know how to hadronize the partons 
on Eout, we can easily add their contribution to dE/dy. 
At ^ = 775 = 0 , we find 

Sout Eont 

= r, [ 


dy 


L 


d^<T^ 


dris 


.(38) 

Fout 


Here f{x,p) is the distribution function for the initial 
free-streaming partons, and we again used the property 
f{x,p) r\j S{y — rjs) to convert 

d^p±d^a^ = {dyd^p±)^^^ = n^{dy d'^pj_){Tsd'^rj_). (39) 

The final total energy after free-streaming and hydrody- 
namical evolution is the sum of the contributions from 
Eout the hydrodynamic freeze-out surface Efo. We 
can now rescale the initial gluon distribution function 
such that for each switching time Tg the final total en¬ 
ergy reproduces the “standard” value defined above. We 
repeat that holding the final energy dE/dy fixed is not 
equivalent to demand ing fix ed final multiplicity dNch / dr] 
(see discussion in Sec. IV A). This should be kept in mind 
when interpreting the results from the following study of 
the sensitivity of physical observables to the switching 
time Tg (i.e. to the time the system needs to thermalize 
sufficiently for hydrodynamics to become applicable). 

Because of our inability to convert the partons on Eout 
to hadrons, we do not know how to include their contri¬ 
bution in final hadronic observables such as transverse 
momentum spectra and flow anisotropies. These observ¬ 
ables are computed by only including hadrons emitted 
from Efo'. 


dNj 

dyd(t)p 


J p±dpi_ J^p^d^ap{x) fi{x,p), 
W r d6 finelyp 

V r d6 ’ 

l^i ^Ypdydct)p 


(40) 

(41) 


where the sum over i runs over all hadron species. Ob¬ 
viously, theoretical predictions based on this procedure 
are not trustworthy if the neglected contribution from 
Eout presents a significant fraction of the total energy. 
Fig-i shows that the associated error increases rapidly 
with Tg, but remains below 5% for 3 fm/c. If the 
system takes longer than 3 fm/c to thermalize, we can 
no longer ignore the “corona” m of particles emerging 
from Eout • We will therefore consider comparisons of our 
predictions for hadronic observables with experimental 
data “meaningful” only for runs with ^3 fm/c. 


found a 9% discrepancy between the value of 11 obtained by re¬ 
constructing it from (i/bulk,i using its kinetic definition [^3 and 
the value obtained directly by applying the projection ( |11| ) on 
the hydrodynamic energy-momentum tensor. Since in our case 
effects from 11 on the hadron spectra and ffow coefficients, once 
integrated over the entire freeze-out surface, are small, we here 
ignored this 9% discrepancy. 
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FIG. 8. (Color online) The fraction (in percent) of the total 
hnal energy ^ contributed by non-thermalized particles on 
Eont, as a function of the switching time Ts- 


VI. ENERGY FLOW ANISOTROPY 


Due to the lack of the contribution from partons on 
Eout, foi* large switching times we are unable to fully re¬ 
construct the hadron azimuthal distribution and its flow 
anisotropy. This is a big handicap when trying to explore 
the effect of large switching times on those physical ob¬ 
servables: The largest effects one sees are due to the loss 
of particles through the Eout surface. 

One way around this problem is to construct the flow 
anisotropy from the azimuthal distribution of energy in¬ 
stead of that of the hadrons. Since the energy angular 
distribution closely tracks that of the particles, the en¬ 
ergy flow anisotropies should be strongly correlated with 
the (momentum-integrated) hadron flow anisotropies. In 
this section, we construct the energy flow anisotropy and 
calibrate it by comparing it to the hadron flow anisotropy 
for small switching times, when the contributions from 
Eout are negligible. 

The azimuthal distributions of energy emitted from 
Eout and Efo at ^ = 77 ^ = 0 are given by 


dE 


dyd(f)p 

dE 


^fo 


EtT^ / n''puP±dpj_ [ p^(fcT^fi, (42) 
“ ( 27 r y Jo Jy, 


dyd(j)r, 


9Ts 

( 27 r)' 


[ n’'p^p]_dp± f d^r±f. 

do t/yiout 


(43) 


In (p2|) i runs over all hadron species, while (43) contains 


only the gluons generated from the KLN model. In ( [43| ) 
we used that on Eout 5 P^d^(7p=p^Tsd?r^ for massless 
partons at ^ = 0. Note that at ^ = 0, rEpj, 
massless partons in (43) while rEp^ = m_L,i = \/ 
for the hadron species i in (42). 



FIG. 9. (Golor online) Scatter plot of the hadron flow 
anisotropies 772,3 for pions (top panels), kaons (middle 
panels) and protons (bottom panels), and energy flow 
anisotropies 1 ^ 2,3 from 400 fluctuating hydrodynamic events 
using Ts=0.6fm/c for the switching time. The left (right) 
column is for elliptic (triangular) flow. The slope and correla¬ 
tion coefficients from a linear ht to the scatter plots are noted 
in each panel. 


Summing these two contributions to a single distribu¬ 
tion dE/dyd(p^ the Fourier coefficients of this azimuthal 
energy distribution can be extracted by the same pro¬ 
cedure as used for calculating hadronic anisotropic flow 
coefficients: 


(44) 


Wn is the energy flow anisotropy coefficient, which quan¬ 
tifies the azimuthal distribution of the energy contributed 
by all the particles - non-thermalized partons as well as 
frozen-out, hydro dynamically flowing hadrons. is the 

energy flow angle associated with Wn- 

We will now show that, for small whose the con¬ 
tributions from Eout to energy and particle emission are 
negligible, Wn from Eq. (44) and Vn from Eq. (41) are 
tightly correlated. The advantage of Wn over Vn is that 
it is easy to include all contributions, including that from 
Eout, in onr calculation of the energy flow, while Vn only 
accounts for contributions from Ej^ and thus misses a 
large fraction of the emitted hadrons when is large. 
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FIG. 10. (Color online) Mean transverse momentum for pions, kaons and protons, with and without pre-equilibrium dynamics, 
compared with data from the ALICE collaboration [32] (horizontal bands), for both MC-KLN (a) and MC-Glauber (b) initial 
models. 


If Wn and Vn are tightly correlated for small where 
both Vn and account for essentially all emitted parti¬ 
cles, we are allowed to use Wn as a proxy for the true Vn 
also for large where Vn computed from (41) no longer 
faithfully represents the full system. 

Fig.i shows the correlation between the hadron flow 
anisotropies vl^ for three selected particle species i and 
the energy flow anisotropies Wn for harmonic orders n = 2 
and 3. For this comparison we chose Ts =0.6 fm/c to 
guarantee that the contribution to Wn from Eout is neg¬ 
ligible. We observe almost perfect correlations, with cor¬ 
relation coefficients very close to 1, for both W 2 vs. V 2 
and ws vs. v^. The slopes w^/v^ are larger than 1 for all 
hadron species i, demonstrating a stronger sensitivity of 
the energy flow coefficients to the hydrodynamic flow 
anisotropies than of v^ for individual hadrons. Among 
the hadrons, heavier species such as protons are more 
sensitive to hydrodynamic flow than lighter species [T], 
but even for protons 'r’2,3 are still smaller than 1^2,3. We 
confirmed similarly strong correlations between ^’2,3 and 
1^2,3 at other switching times < 2.5 fm/c but saw that 
the correlation gradually breaks down for large values 
Tq >3 fm/c when too much of the total energy emerges 
from Eout- 


VII. CONSTRAINING THE DURATION OF 
THE PRE-EQUILIBRIUM STAGE 

In this section, we will use the energy flow anisotropies 
rc2,3 and the mean transverse momenta for pions, kaons 
and protons to constrain the duration of the pre¬ 


equilibrium stage. It is well known that the average 
transverse momenta of hadrons with different masses 
help to separate random thermal motion (i.e. the freeze- 
out temperature) from the effect of collective radial hy¬ 
drodynamic flow in the final state. Switching from initial 
free-streaming to hydrodynamics at later times means 
more initial flow after Landau matching but less time for 
developing hydrodynamic flow. Fig. has already shown 
that the net effect is an increased radial flow at freeze-out 
which will lead to harder momentum spectra and an in¬ 
crease in the average px • If is therefore expected that the 
measured (p^) values for pions, kaons and protons will 
put an upper limit on the switching time, by limiting the 
amount of radial flow at freeze-out. 

This is illustrated in Fig. [Tq] where the mean transverse 
momenta (p^) for pions, kaons and protons are plotted 
as a function of switching time r^, with MC-KLN ini¬ 
tial conditions (propagated with specific shear viscosity 
r]/s = 0.2) in panel (a) and MC-Glauber initial conditions 
(propagated with 77/5 = 0.08) in panel (b). Shown for 
comparison are experimental data from the ALICE Col¬ 
laboration [32] that were obtained by extrapolating the 
measured spectra to the full p± range before calculating 
the mean. When pre-equilibrium dynamics is included 
in the calculations (solid lines with filled symbols), the 
mean transverse momenta are seen to increase with r^, 
as anticipated. The effect is strongest for protons whose 
large mass makes them most susceptible to flow. One 
sees that, with the chosen values for the shear viscosity 
and freeze-out temperatures, the MC-KLN initializations 
with free-streaming pre-equilibrium dynamics have diffi¬ 
culties accommodating the data unless one postulates es- 
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sentially instantaneous thermalization, and even then the 
pion mean p±_ is still statistically significantly too large 
(by about 10%) compared to the measurements. On the 
other hand, the smaller shear viscosity used for evolving 
the MC-Glauber initial profiles in panel (b) reduces the 
transverse shear stress and thus builds less radial flow, 
giving room for some pre-equilibrium radial flow. With 
this combination of initial conditions and shear viscosity, 
a switching time around Ifm/c appears to be preferred 
over significantly smaller and larger Tq values. 

In Fig. we also show for comparison as dashed lines 
with open symbols the corresponding results for hydro- 
dynamic evolution without pre-equilibrium dynamics. In 
this case Tq has the meaning of the starting time for the 
hydrodynamic evolution, but with initial conditions that 
have not evolved in the transverse plane between tq and 
Tg. Without pre-equilibrium, delaying the start of the 
hydrodynamic expansion leads to a reduction of the final 
radial flow, since less time is available for its generation 
before the matter reaches the decoupling temperature. 
As a result, (p^) decreases with increasing r^, and the 
effect is again stronger for protons than for pious and 
kaons, due to their larger mass. 

For the purely hydrodynamic runs without pre- 
equilibrium dynamics, we see in Fig. [^that no choice of 
Ts can reproduce all three measured values simultane¬ 
ously, for either of the two initial condition models and 
associated shear viscosities. The reader may wonder how 
this can be consistent with successful earlier fits of the 
measured transverse momentum spectra for these three 
particle species [33l |34]. Part of the answer is that in 
[33l l34] the quality of the model description of the data 



FIG. 11. (Color online) Elliptic and triangular energy flow 
anisotropy coefficients W 2 ^ W 3 as a function of the switch¬ 
ing time from single-shot hydrodynamic simulations with a 
smooth ensemble-averaged MC-KLN initial profile. 


was judged by the overall shape of the p_L-distributions 
whereas here we compare with only a single moment of 
that distribution which, however, measured with very 
good precision. We will return to this question in the 
following section where we look at a somewhat larger set 
of experimental data that are, in addition to the mean 
p±, also sensitive to the shape of the p^-distributions and 
try to fit them by simultaneously varying several hydro- 
dynamic parameters. 

For anisotropic flow the dependence is more sub¬ 
tle: as shown in Fig. a free-streaming pre-equilibrium 
stage reduces the source eccentricity at the start of the 
hydrodynamic evolution, so it reduces the amount of flow 
anisotropy that can be generated during the hydrody¬ 
namic stage in response to this initial eccentricity. How¬ 
ever, as shown in Fig. I§i, it also creates a non-zero hy¬ 
drodynamic flow anisotropy at the start of the fluid stage 
which gives the hydrodynamic evolution of anisotropic 
flow a boost. Fig. pT] shows that, for switching times up 
to about 2 fm/c, the combined effect are final elliptic and 
triangular flow anisotropies W 2,3 that are almost inde¬ 
pendent of the switching time. In fact, both W 2 and W 3 
slightly increase with increasing switching time until Tq 
reaches about 1.5 fm/c. Only for larger switching times 
beyond 2 fm/c does the reduction of ^2,3 before the start 
of the hydrodynamic evolution cut into the finally es¬ 
tablished anisotropic flow coefficients, and for very large 
switching times both W 2 and W 3 approach zero. 

The naive expectation that the pre-equilibrium dilu¬ 
tion of the source eccentricity before Tq should reduce the 
finally established anisotropic flow mils] is borne out 
only if one completely ignores the position-momentum 
correlations created by the pre-equilibrium dynamics and 
starts the hydrodynamic stage with zero transverse flow. 
This is illustrated by the dashed lines with open symbols 
in Fig. 

Figure thus leads, in disagreement with the ear¬ 
lier statements made in Refs. D US], to the (revised) 
conclusion that elliptic and higher-order anisotropic flow 
measurements alone cannot put a tight upper limit on 
the duration of the pre-equilibrium stage. Fig. sug¬ 
gests that, as long as pre-equilibrium contributions to the 
final flow pattern are consistently accounted for, the fi¬ 
nal anisotropic flows are insensitive to how strongly the 
medium is coupled during the first 2 fm/c or so. Whether 
this stage is described hydro dynamically (very strong 
coupling) or by free-streaming partons (very weak cou¬ 
pling) , one sees the same final flow anisotropy. This find¬ 
ing supports the idea of “universal transverse flow” dur¬ 
ing the earliest stages of the fireball evolution proposed 
by Pratt and Vredevoogd [35]. In contrast, radial flow 
(which affects the slope and mean p^ of the transverse 
momentum spectra) exhibits a strong and monotonic Tq 
dependence already for small switching times that can be 
used much more effectively to put an upper limit on the 
thermalization time in heavy-ion collisions. 
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FIG. 12. (Color online) Parameter search result for the KLN 
initial-state model, evolved in single-shot hydrodynamic mode 
following a free-streaming stage of duration Ts—tq. (a) His¬ 
togram of the distribution, (b-d) 2-dimensional projec¬ 
tions of those parameter triplets (ts, 77 /s, Tdec) corresponding 
to < 50. The size of the rings around their positions in¬ 
creases with decreasing x^, be. with increasing fit quality. 


VIII. PARAMETER OPTIMIZATION 


In addition to the switching (or starting) time our 
hydro dynamical model has several other input param¬ 
eters whose choice influences the final physical observ¬ 
ables. In the preceding sections we only varied leav¬ 
ing these other parameters unchanged, in order to gain 
generic insights into which of the different observables 
at our disposal provide the strongest constraints on r^. 
However, it is immediately obvious that there should be 
some sort of tradeoff between effectively weakening the 
interactions in the pre-equilibrium stage, say by length¬ 
ening the free-streaming period, and weakening the in¬ 
teractions during the later hydrodynamic stage, say by 
shortening the free-streaming stage and increasing the 
shear viscosity during the subsequent hydrodynamic evo¬ 
lution. The effects of changing the transition time be¬ 
tween free-streaming and hydrodynaming evolution and 
of changing the shear viscosity during the hydrodynamic 
evolution are therefore entangled, and we should optimize 
both parameters simultaneously. Furthermore, since the 
slopes of the final spectra are controlled by a combina¬ 
tion of the temperature and radial flow on the freeze- 
out hypersurface, and a change in viscosity affects the 
transverse pressure gradients and thus the radial flow, 
we should allow T^ec to vary together with rj/s. 

Here we report on a study where we allow r^, 77/s and 
Tdec to vary simultaneously, trying to find the best com¬ 
bination by comparing the model predictions for 
{p±) 7 r+i {p±)k+i and {pi_)p for 2.76 ATeV Pb+Pb colli¬ 
sions of 10-20% centrality with experimental data from 




FIG. 13. Gomparison with experimental data of single-shot 
hydrodynamic model predictions with the best ht parame¬ 
ter combination (smallest x^ value) for KLN-initialized sim¬ 
ulations with pre-equilibrium free-streaming dynamics, for 
2.76 A TeV Pb+Pb collisions at 10%-20% centrality. The 
best-fit parameter values are listed in panel (a), (a) Trans¬ 
verse momentum spectra for tt^, and protons, compared 
with ALIGE data [JO]. (b) Eccentricity-scaled -differential 
elliptic and triangular flow anisotropies for charged hadrons, 
compared with ATLAS data m- Single-shot hydrodynamic 
simulations predict the ensemble averaged flow Vn which was 
scaled by the corresponding ensemble-averaged eccentricity 
Cn (solid lines). The ATLAS anisotropic flow data are mea¬ 
sured with the event-plane method and are thus affected by 
the variance of the event-by-event Vn distribution. They were 
fitted with a smooth curve and scaled by the rms eccentric¬ 
ity €n{2} (dashed lines). The shaded area around the dashed 
lines represents the experimental error of the ATLAS UnjEP} 
measurements m- 


the ALICE [32] and ATLAS |36| collaborations and min¬ 
imizing the value of x^- The values of the experimental 
measurements for these observables are summarized in 
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FIG. 14. (Color online) Same as Fig. 12 but for single-shot 
hydrodynamic simulations starting at Ts without preceding 
pre-equilibrium stage. 


Table H 




0.0782 

± 

0.0019 



0.0316 

± 

0.0008 

(px)^ 

+ (GeV/c) 

0.517 

± 

0.017 

{ p ±) k + (GeV/c) 

0.871 

± 

0.027 

(Px)p 

(GeV/c) 

1.311 

± 

0.034 


TABLE I. Experimental data for the five hadronic observables 
from 2.76 A TeV Pb+Pb collisions of 10-20% centrality that 
were considered in our fit. The mean -integrated elliptic 
and triangular flow values for charged hadrons were measured 
by ATLAS [36], the mean transverse momenta for positively 
charged pions, kaons and protons (extrapolated to the full p± 
range) by the ALICE Collaboration [3^ . 


A somewhat more ambitious fit with five hydrody¬ 
namic model parameters and three experimental observ¬ 
ables (charged multiplicity, (^2^) and {vf^)) at six col¬ 
lision centralities each (i.e. alltogether 18 observables), 
including a hadronic afterburner but no pre-equilibrium 
dynamics, was recently reported in [37]. We perform 
simulations both with and without pre-equilibrium dy¬ 
namics, in order to assess its impact on the best-fit val¬ 
ues for the other model parameters. Our simulations are 
done in single-shot mode with smooth ensemble-averaged 
initial density profiles, not in event-by-event mode with 
fluctuating initial profiles as the work reported in 
Since we used data on the mean elliptic and triangular 
flows of charged hadrons, obtained by ATLAS [36] from 
their full reconstructed event-by-event probability distri¬ 
butions, instead of their rms values that were used in 
and which are affected by the variance of their event-by¬ 
event fluctuations, the added numerical cost of event-by¬ 
event hydrodynamic simulations could be avoided in our 




EIG. 15. Same as Eig. [^ but for single-shot hydrodynamic 
simulations starting at Ts without preceding pre-equilibrium 
stage. 


analysis. 

The parameter space is explored by Latin hypercube 
sampling [38] |39] , which is a statistical sampling method 
for optimizing the selection of parameter combinations 
from a high-dimensional parameter space. Because our 
parameter space is 3-dimensional, each point drawn from 
this space by the Latin hypercube method contains three 
components. In order to efficiently span the parameter 
space, no two points in the sample share the same value 
for any of the three parameters. We checked that, af¬ 
ter projecting the sampled triplets onto any of the three 
components, the distribution of that component was uni¬ 
form. 

After running the hydrodynamic simulation for a given 
triplet of parameters from the Latin hypercube sample, 
the quality of the description obtained with this parame¬ 
ter set is estimated by computing the of fho resulting 
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FIG. 16. (Color online) Same as Fig. 12 but for MC-Glauber 
initial conditions, free-streamed to Ts and then evolved hy- 
drodynamically. 


fit of the selected experimental data: 


(45) 


Here the sum runs over the five observables, Ei is the 
value and the combined statistical and systematic 
error of the experimental measurement of the observ¬ 
able, and Oi(Ts, 77/s, Tdec) is the value of the observable 
from the simulation with parameter set (r^,77/s,Tdec)- 
Since we are using single-shot hydrodynamics with an 
ensemble-averaged initial profile, Oi has no statistical er¬ 
ror. Assuming that the five chosen observables are un¬ 
correlated, with three fit parameters we have two sta¬ 
tistical degrees of freedom, and a good fit of the data 
should thus have x^/2:^l. The best fits we have been 
able to achieve with the three selected model parameters 
have x^/2 ^ 9 — 15 (see Table III). This suggests that 
at least one additional physical mechanism not captured 
by these three model parameters may play an important 
role in describing the chosen observables. This could be, 
for example, that hydrodynamics breaks down during the 
final hadronic stage of the fireball expansion and needs 
to be replaced there by a microscopic model for hadronic 
rescattering nain]. We have not tested this hypothesis. 

For the KLN initial-state model we drew a Latin 
hypercube sample of 1300 points covering the fol¬ 
lowing parameter ranges: 0.1 fm/c < < 1.4 fm/c, 

0.08 < 77/5 < 0.28, and 100 MeV < T^ec < 170 MeV. For 
the Glauber initial-state model we sampled 950 points 
in the range 0.1 fm/c < Tg < 2 fm/c, 0 < 77 /g < 0T6, and 
100 MeV < Tdec < 170 MeV0Figs. [l2) [HI [l6]and[l8jshow 



FIG. 17. Same as Fig. but for MG-Glauber initial condi¬ 
tions, free-streamed to Ts and then evolved hydrodynamically. 


the corresponding distributions for simulations with 
KLN and Glauber initial conditions with and without 
a free-streaming pre-equilibrium stage, respectively (see 
captions). The best-fit values for the three model pa¬ 
rameters in each of the four cases, and the predictions 
the model makes for these best-fit parameter values for 
the observables listed in Table [Tl are summarized in Ta- 


show scatter plots of all simulations with a total < 50 
in each of the three 2-dimensional projections of the 3- 
dimensional parameter space. The size of each blob indi¬ 
cates the quality of the description of the data achieved 
with the corresponding parameter set: The larger the 
blob, the smaller the total x^ and the better the fit. 


bles \n\ and |III[ Panels (b-c) in Figs. 


12 


14 


16 


and 


18 


2 In both cases the upper end of the explored range for Ts is small 


enough that the particle and energy losses through the corona 
studied in Sec. ^ can be ignored. 
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FIG. 18. (Color online) Same as Fig. 16 but for single-shot 
hydrodynamic simulations starting at Ts without preceding 
pre-equilibrium stage. 


Figures and show the distributions of 

pious, kaons and protons (a) and of the charged hadron 
elliptic and triangular flows predicted by the simulations 
with the best-fit parameter sets, for each of the four sim¬ 
ulation modes listed in Tables [TTj m We use the equa¬ 
tion of state s95p-PCE EH El which assumes chemi¬ 
cal freeze-out of hadron abundances at Tchem = 165 MeV 
and therefore overpredicts the measured proton yields 
by about 50% (without the hadronic afterburner, our 
hydrodynamic approach cannot account for baryon-anti- 
baryon annihilation during the final hadronic rescattering 



Model 

Pre-eq. 

Ts (fm/c) 

v/s 

Tdec (MeV) 

MG-KLN 

Yes 

0.129 

0.206 

108 

MG-KLN 

No 

0.131 

0.186 

113 

MG-Glb 

Yes 

0.597 

0.160 

111 

MG-Glb 

No 

0.295 

0.147 

104 


TABLE II. Best-fit parameters from a parameter search for 
four different types of simulations. 


Model MC-KLN MC-KLN MC-Glb MC-Glb 

Pre-eq. Yes No Yes No 



0.083 

0.088 

0.070 

0.071 


0.030 

0.030 

0.034 

0.034 

{p±}^+ (GeV/c) 

0.550 

0.545 

0.539 

0.518 

{P±}k+ (GeV/c) 

0.900 

0.877 

0.869 

0.852 

{p±}p (GeV/c) 

1.349 

1.302 

1.293 

1.279 


18.624 

33.590 

29.541 

23.931 


TABLE III. Model predictions for the five observables listed 
in Table 1^ with the best-fit parameters listed in Table |TT] 


FIG. 19. Same as Fig. but for single-shot hydrodynamic 
simulations starting at Ts without preceding pre-equilibrium 
stage. 


stage, which is required to reproduce the experimental 
yields [l5]). For the protons one should therefore ignore 
the normalization of the p± spectra and focus instead on 
their shape. Clearly, the description of the experimental 
data in Figs. 13, 15, 17 and (of which only the lowest 
non-trivial moments were used in the fit) is not perfect, 
but of similar quality as most other, less systematic pa¬ 
rameter fits published in the literature. 


At first sight, though, there appears to be one ex¬ 
ception: As seen in Fig. [T9] 3, the MC-Glauber model 
without pre-equilibrium dynamics appears to provide a 
qualitatively better simultaneous description of the dif¬ 
ferential charged hadron elliptic and triangular flows than 
the other approaches. (In particular, the KLN model 
has troubles to describe these two observables simulta¬ 
neously, as has been noted before [46].) However, the 
overall of “Glauber without pre-equilibrium” fit 
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is worse than that for the “KLN with pre-equilibrium” 
fit, due to a larger discrepancy with the data of the pre¬ 
dicted mean p± for kaons and pious. Fig. [T^ also shows 
that the best “Glauber without pre-equilibrium” simu¬ 
lation does not describe the slope of the pion spectrum 
quite as well as the best “KLN with pre-equilibrium” fit 
in Fig. [T^. This clearly detracts from the apparently 
much better description of the px-differential V 2 %{p±) in 
Fig. [T^ compared to Fig. [T^ which may simply reflect 
an incorrect weighting of regions of high and low p± in the 
best-fit “Glauber without pre-equilibrium” simulation. 

Let us discuss a few other trends that are visible in 
Figs T 2 ][T^ First, panels b and c in Figs. \TE 

and m as well as Table |TTj show that the inclusion of 
pre-equilibrium dynamics puts some upward pressure on 
the best-fit value of 77 /s, for both types of initial condi¬ 
tions. The effect is not large (of order 10%) but appears 
to be significant. (We will not be absolutely certain of 
the significance of this observation until we have com¬ 
pleted a full Markov Chain Monte Carlo (MCMC) sim¬ 
ulation of the posterior model parameter distributions 
[371117] which is beyond the scope of this paper.) Sec¬ 
ond, by looking at the spread of the blobs shown in panel 
b of Figs. and we see a reduced sensitiv¬ 

ity of the fit quality to the decoupling temperature when 
pre-equilibrium dynamics is accounted for in the simula¬ 
tions: the same upper limit of 50 for the total allows 
for larger deviations of Tdec from the best-fit value if the 
simulations include a pre-equilibrium stage. Third, it can 
be seen from panel c in the same figures that allowing for 
pre-equilibrium build-up of flow reduces the probability 
of a good fit for switching times that significantly ex¬ 
ceed the best-fit value in the KLN model but not in the 
Glauber model. This is consistent with Fig. where 
we observed that allowing for any appreciable delay of 
the hydrodynamic stage due to pre-equilibrium dynam¬ 
ics tends to cause the model to overpredict the mean p± 
of the hadrons (especially the protons) for KLN initial 
profiles (panel a) but not for the Glauber model profiles 
(panel b). Indeed, for the Glauber model the neglect of 
pre-equilibrium flow puts downward pressure in r^, be¬ 
cause of the effect of missing radial flow on the mean 
hadron p±. One sees this in both Fig. and when 
comparing panels c in Figs. [Tb] and [Tsl 


IX. SUMMARY AND CONCLUSIONS 

In this paper, we thoroughly investigated the effect of 
pre-equilibrium flow on heavy-ion collision observables. 
In order to estimate the maximum difference between 
simulations where the hydrodynamic stage was initiated 
with non-zero flow resulting from a preceding stage of 
pre-equilibrium evolution and others where all dynamical 
effects before the onset of hydrodynamic expansion were 
simply ignored (as has been common practice in many 
earlier studies), we here assumed the pre-equilibrium 
stage to be non-interacting, i.e. free-streaming. By ad¬ 


justing the switching time between the free-streaming 
and hydrodynamic stages we can thus smoothly switch 
between a picture where the earliest stage of the collision 
is coupled infinitely weakly and one were it is coupled 
infinitely strongly. 

The pre-equilibrium and hydrodynamic stages are 
Landau-matched, conserving energy and momentum and 
fully accounting for all components of the pre-equilibrium 
energy-momentum tensor with non-zero (and possibly 
large) dissipative pressure components (bulk and shear 
viscous pressures). The study presented here was done 
for zero bulk viscosity, i.e. the non-zero bulk viscous 
pressure generated by the Landau matching procedure 
was allowed to evolve dynamically to zero with a short 
microscopic relaxation time. The viscous shear stress at 
the beginning of the hydrodynmic stage was found to 
be large, corresponding to an inverse Reynolds number 
slightly above unity for early switching times and further 
increasing if one switches from free-streaming to hydro¬ 
dynamics later. In the hydrodynamic stage, these large 
starting values for the shear stress decrease quickly on 
a short microscopic time scale, due to the assumed low 
shear viscosity 77 /s of the hydrodynamic fluid. 

On the way we had to solve the problem that extended 
free-streaming leads to a loss of a significant fraction of 
the energy of the expanding fireball because the volume 
of the “corona” of dilute matter that never thermalizes 
and therefore never becomes part of the hydrodynamic 
fluid increases. While we are not able to convert corona 
partons into hadrons we can account for their energy, so 
we can include them in anisotropic flow coefficients 
that characterize the azimuthal anisotropy of the energy 
flow even if we cannot account for their contributions 
to anisotropic particle flow Vn. We demonstrated that 
for small switching times were corona losses can be ne¬ 
glected, Wn and Vn closely track each other. 

Our most important finding is that, contrary to tradi¬ 
tional believes (previously also held by one of the present 
authors), the anisotropic flow coefficients are not very 
sensitive to the thermalization (or switching) time as 
long as the latter does not significantly exceeed about 
2 fm/c. The reason is that anisotropic flow not generated 
in the hydrodynamic stage, due a reduced initial fireball 
eccentricity and shortened hydrodynamic stage when 
is allowed to grow large, is almost precisely compensated 
for by anisotropies in the space-momentum correlations 
that develop during the free-streaming phase and man¬ 
ifest themselves as non-zero starting values for the mo¬ 
mentum anisotropy of the energy-momentum tensor in 
the hydrodynamic stage. On the other hand, we found 
that an extended weakly coupled pre-equilibrium stage 
leads to stronger radial flow at the end of the hydrody¬ 
namic evolution, caused by large initial flow at the begin¬ 
ning of the hydrodynamic stage which is stronger than 
it would have been if the pre-equilibrium stage had been 
strongly coupled (i.e. describable by hydrodynamics). 
This is a result of the faster signal propagation speed 
in the pre-equilibrium stage which, in the free-streaming 
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case, is given by the speed of light and thus exceeds the 
hydrodynamic speed of sound by almost a factor 2. This 
kick-start of the radial flow by pre-equilibrium evolution 
overcompensates the loss of radial flow generated dur¬ 
ing the hydrodynamic stage which is shortened when the 
pre-equilibrium stage lasts longer. 

It thus turns out that the mean transverse momentum 
{p±) of the finally emitted hadrons, which (especially for 
the massive baryons) is strongly affected by the final ra¬ 
dial flow of the fireball, provides a stronger upper limit on 
the thermalization time than the anisotropic flow coef¬ 
ficients. This agrees with recent findings by Romatschke 
and collaborators [6]. 

Of course, the switching time is only one of several 
model parameters affecting the final observables. In the 
last section we therefore performed an extended param¬ 
eter search where we varied the switching time, decou¬ 
pling temperature and specific shear viscosity of the fluid 
simultaneously, for two different initial state models and 
both runs with and without pre-equilibrium dynamics. 
While this does not exhaust the list of parameters and 
possibilities this exercise provides a more holistic picture 
of the effects of pre-equiibrium dynamics on final observ¬ 
ables and on the values of medium parameters extracted 
from a comparison of the model predictions with exper¬ 
imental data. We found for both initial state models 


that accounting for pre-equilibrium dynamics slightly in¬ 
creases the optimal values for the specific shear viscosity 
r]/s extracted from mean p± and anisotropic flow mea¬ 
surements and reduces the sensitivity of the observables 
to the decoupling temperature. However, the extracted 
limits for the thermalization time turned out to depend 
sensitively on the model for the initial density profiles. 
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